• Nested ANOVA or nested design
  • factors might be hierarchical - in other words nested - within one another
  • The sources of variance are therefore hierarchical too
• The factorial ANOVA design is the most common experimental design used to investigate more than one treatment variable
  • In a factorial design every combination of treatments from two (or more) treatment variables is investigated.
  • The main purpose of a factorial design is to evaluate possible interactions between variables.
  • An interaction between two explanatory variables means that the effect of one variable on the response depends on the state of a second variable.

• Nested designs are hierarchical
  • often contain sub-replicates that are random, uncontrolled, nuisance effects
  • but the nested factors can be of interest too
• Factorial designs are
  • all pairwise combinations,
  • and often involve all combinations of factor levels
  • when each factor is fixed interactions can be assessed
• Completely nested designs therefore have no interaction terms, whereas factorial designs do
• Mixed models can have a combination of fixed and random factors that are more complicated

• Example 1: Study of “repeatability” (simple nested design)
• The walking stick, Timema cristinae, is a wingless herbivorous insect on plants in chaparral habitats of California.
• Nosil and Crespi (2006) measured individuals using digital photographs.
• To evaluate measurement repeatability they took two separate photographs of each specimen.
• After measuring traits on one set of photographs, they repeated the measurements on the second set.

andrew_data <- read.table('andrew.tsv', header=T, sep=‘\t')
head(andrew_data)

• There are four variables: ‘TREAT’, ‘PATCH’, ‘QUAD’ and ‘ALGAE’
• The main effect factor is TREAT
• Make a simplified factor called TREAT2, in which 0% and 33% are a level called “low” and 66% and 100% are “high”

andrew_data$TREAT2 <- factor(c(rep(“low”,40),rep(“high”,40))

• The nested factor is PATCH - also need to turn this into a factor

andrew_data$PATCH <- factor(andrew_data$PATCH)

• In this case, our response variable is ALGAE
• Look at the distribution of ALGAE for the two levels of TREAT2 using boxplots based on the patch means, which are the replicates in this case.

andrew.agg <- with(andrew_data, aggregate(data.frame(ALGAE), 
                  by = list(TREAT2=TREAT2, PATCH=PATCH), mean)

library(nlme)
andrew.agg <- gsummary(andrew_data, groups=andrew_data$PATCH)

boxplot(ALGAE ~ TREAT2, andrew.agg)

• Evaluate assumptions based on the boxplots
• Is the design balanced (equal numbers of sub-replicates per PATCH)?

• Run the nested ANOVA:

nested.aov <- aov(ALGAE ~ TREAT2 + Error(PATCH), data=andrew_data)
summary(nested.aov)

• Do we detect an effect of TREAT2 (high vs low sea urchin density)?
  
• Estimate variance components to assess relative contributions of the random factors

library(nlme)
VarCorr(lme(ALGAE ~ 1, random = ~1 | TREAT2/PATCH, andrew_data))

• Calculate the % of variation due to between-treatment differences vs. due to among patches within treatment differences.
• See pg. 302 in Logan if you need help.
• What do these variance component estimates tell us???

• For example, Relyae (2003) looked at how a moderate dose (1.6mg/L) of a commonly used pesticide, carbaryl (Sevin), affected bullfrog tadpole survival.
• In particular, the experiment asked how the effect of carbaryl depended on whether a native predator, the red-spotted newt, was also present.
• The newt was caged and could cause no direct harm, but it emitted visual and chemical cues to other tadpoles
• The experiment was carried out in 10-L tubs (experimental units), each containing 10 tadpoles.
• The four combinations of pesticide treatment (carbaryl vs. water only) and predator treatment (present or absent) were randomly assigned to tubs.
• The results showed that survival was high except when pesticide was applied together with the predator.
• Thus, the two treatments, predation and pesticide, seem to have interacted.

continuous response and two main effect variables

rnadata <- read.table('RNAseq.tsv', header=T, sep='')
gene <- rnadata$Gene80
microbiota <- rnadata$Microbiota
genotype <- rnadata$Genotype

make new “pseudo factor,” combining genotype and microbiota

gxm <- interaction(genotype,microbiota)
levels(gxm)
boxplot(gene ~ gxm)

specify the following 2 contrasts

contrasts(gxm) <- cbind(c(2, -1, 0, -1), c(-1, -1, 3, -1))

Fit the factorial linear model

rna_aov <- aov(gene ~ gxm)

Examine the ANOVA table, using supplied contrasts. Figure out the appropriate titles to give them.

summary(rna_aov, split = list(gxm = list('xxx'=1,'xxx'=2)))

What does the contrast summary tell you about the nature of the interaction?

rnadata <- read.table('RNAseq.tsv', header=T, sep='')
head(rnadata)

variables excluding first 5 and last 5 observations

gene <- rnadata$Gene80[6:75] 
microbiota <- rnadata$Microbiota[6:75]
genotype <- rnadata$Genotype[6:75]
boxplot(gene ~ microbiota)
boxplot(gene ~ genotype)
boxplot(gene ~ microbiota*genotype)

Estimate the variance components using Restricted Maximum Likelihood (REML)

library(lme4)
lmer(gene ~ 1 + (1 | microbiota) + (1 | genotype) + (1 | microbiota:genotype))

Based on the REML sd estimates, what are the relative contributions of the factors to total variance in gene expression?

longevity_data <- read.table('longevity.csv', header=T, sep=',')
head(longevity_data)

Variables

long <- longevity_data$LONGEVITY
treat <- longevity_data$TREATMENT
thorax <- longevity_data$THORAX

• check to see if the covariate should be included

boxplot(long ~ treat)
plot(long ~ thorax)

• assess assumptions of normality and homogeneity of variance

plot(aov(long ~ thorax + treat ), which = 1)

• †ry it again with a transformed response variable

plot(aov(log10(long) ~ thorax + treat ), which = 1)

• visually assess linearity, homogenetiy of slopes and covariate range equality

library(lattice)
print(xyplot(log10(long) ~ thorax | treat, type = c("r", "p")))

• formally test homogenetiy of slopes by testing the interaction term

anova(aov(log10(long) ~ thorax*treat))

• formally test covariate range disparity by modeling the effect of the treatments on the covariate

anova(aov(thorax ~ treat))

• FINALLY, set up contrasts, fit the additive model and visualize the results (pg. 459 and 460 of your Logan book)
• Summarize the trends in a nice plot (pg. 461 of your Logan book)